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^ Abstract 

PH The aim of this paper is to prove results about the existence and stabil- 

\) ity of multiple steady states in a system of ordinary differential equations 

_! introduced by R. Lev Bar-Or [J to model the interactions between T cells 

. ^ and macrophages. Previous results showed that for certain values of the 

r^ parameters these equations have three stationary solutions, two of which 

JL^ are stable. Here it is shown that there are values of the parameters for 

1 I which the number of stationary solutions is at least seven and the num- 

ber of stable stationary solutions at least four. This requires approaches 
'"H different to those used in existing work on this subject. In addition, a 

^ rather explicit characterization is obtained of regions of parameter space 

^•-^ for which the system has a given number of stationary solutions. 

^ 1 Introduction 

^—^ Many phenomena in biology and medicine can be modelled by systems of ordi- 

^^ nary differential equations. Sometimes models are used which have very large 

. . numbers of unknowns and parameters. Often quantitative information about 

^ ^ these parameters is obtained from high throughput experimental techniques. 

K> This kind of approach plays a central role in the field of systems biology. (For 

^ an introduction to systems biology see [3].) A complementary approach is to try 

to capture important effects with a description in terms of dynamical systems 
with a few unknowns. At the same time large numbers of parameters may be 
handled by looking for qualitative features of the solutions which are present for 
large regions of parameter space. This leads to mathematical problems more 
accessible to analytical treatment and less dependent on large-scale simulations. 
Some thoughtful remarks on the relative advantages of these different paths to 
understanding in the applications of mathematics to biology can be found in 

One type of qualitative behaviour which is frequently of interest for biologists 
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is that of multistability, i.e. the existence of several stable steady state solutions 
for given values of the parameters. This means that a biological system can 
function in more than one way, depending on its history. In addition to the 
intrinsic interest of this fact it opens up the possibility of manipulating the 
biological system in some way which is favourable for practical applications 
such as therapies for certain diseases. 

The subject of what follows is a system of four ordinary differential equa- 
tions with eighteen parameters introduced in [4] to model certain aspects of 
the immune system. A mathematically rigorous analysis of some featues of the 
behaviour of solutions of this system was carried out in [7]. In particular it 
was proved that for some values of the parameters there is only one stationary 
solution while for other parameter values there are three stationary solutions, 
two of which are stable. When there is bistability the two stationary solutions 
can be identified with two different states of the immune system which are said 
to be Thl-dominated and Th2-dominated, respectively. They are distinguished 
by the concentrations of certain substances called cytokines which immune cells 
use to communicate with each other. In the derivation of the model the cells 
which produce and react to the cytokines are T cells and macrophages. The 
unknowns in the ODE system are concentrations of cytokines and the cells do 
not occur explicitly in the model. 

The function of the immune system is to defend the host against harmful in- 
fluences such as pathogens and toxins which can cause diseases. (For a detailed 
introduction to immunology see [S] or [H].) Under some circumstances the im- 
mune system may malfunction by attacking host tissues, leading to autoimmune 
diseases. The choice between Thl- and Th2-dominated immune responses may 
have an important influence on the course of diseases, whether they are due to 
pathogens or autoimmune in nature. More information on this can be found in 
and [Ij. A survey of work on the development of Thl and Th2 responses in 
an individual organism and approaches to understand this phenomenon using 
mathematical modelling is given in [1]. 

In what follows it is shown that for suitable values of the parameters the sys- 
tem of ^4, has seven stationary solutions, four of which are stable. This means for 
instance that it is possible to have two steady states which are Thl-dominated 
with the degree of dominance being different. This property of the system had 
not previously been observed. As noted in [7], the property of bistability can be 
obtained if only T cells are taken into account. On the other hand the influence 
of the macrophages is essential in order to fulfil the assumptions of the theorems 
in what follows which assert the existence of more than two stable steady states. 

The essential qualitative features of the dynamics exhibited in what follows 
can already be found in a much simpler system which arises in a special case. 
This model system has only two unknowns and two parameters. The basic sys- 
tem of [1] and the two-dimensional model system are introduced in Sect. l2J Sect. 
[3] contains statements and proofs of theorems on multistability for the model 
system. This allows basic ideas to be presented in a relatively simple context. 
Some of the techniques used extend in a straighforward way to the analysis 
of the four-dimensional system. In fact they suffice to prove the existence of 



four stable stationary solutions for an open set of the parameter space which is 
described rather explicitly. To prove the existence of more steady states of the 
four-dimensional system a significant refinement of the techniques is necessary. 
These ideas are presented in Sect. |4j The last section is devoted to further 
discussion of the results of the paper and possible generalizations. 

2 The basic equations 

The basic dynamical system studied in this paper is 

dx 

-^ = -d,Xi + g(hi); i = 1,2,3,4. (1) 

The di are positive constants. The function g is given by 

5(x) = i(l+tanh(x-0)) (2) 

where ^ is a constant. The functions hi are defined by hi = '^- aijXj for some 
constants aij. The equations were written in this form in [7]. The system 
was first introduced in [3] in a different notation. The relation between the 
two notations is explained in [7j. The coefficients in (IT]) are required to satisfy 
a number of conditions which will now be listed. The formulation of these 
conditions has been changed slightly from that used in [7] so as to make certain 
calculations more efficient. In the formulation used here each aij is of the form 
hij + Cij and the following conditions hold: 

1. {-if+'h, > 

2. bij — —b2j and 63^ = —b^j for all j 

3. c^j — C4j = for all j 

4. there is a constant A' > such C2j = Kcij for all j. 

5. cii > and C13 > 

The last of these conditions is a slight weakening compared to [4] , where strict 
inequality was required. The analysis done in what follows only requires the 
assumptions 1. and 3. Conditions 2., 4. and 5. have been included here only to 
make the relation to the original set-up of ^ clear. 

The interpretation of these equations is as follows. The quantities xi and X3 
are the concentrations of Thl cytokines produced by T cells and macrophages 
respectively. The quantities X2 and x^ are the corresponding concentrations of 
Th2 cytokines. Thus zi — xi-\- x^ and Z2 = X2 + x^ are the total concentrations 
of Thl and Th2 cytokines, respectively. The signs of the bij reflect the fact 
that each type of cytokine (Thl, Th2) stimulates cells to produce that type 
of cytokine and inhibits their production of the other type. The coefficients 
Cij encode the effect of antigen presentation. The cytokines modify the way 



in which macrophages present antigen to T cehs and this in turn modifies the 
cytokine production of the T cells. The signs of these coefficients are based on 
experimental data reviewed in |4j . Setting the coefficients Cij to zero corresponds 
to neglecting the effect of antigen presentation. 

A special case of this system which is important in what follows is obtained 
by requiring that di = 1 for all i, 9 — 0, the quantities \bij\ with j odd are all 
equal, the quantities \bij\ with j even are all equal and the quantities Cij with 
i = 1,2 are all equal. Then defining A = bu, B = 612 and C = cn gives a 
system which is (21)-(24) of JJ. It implies a closed system for the variables zi 
and Z2 which is (25)-(26) of [7]. The special case A ^ B oi the last system is 
given by 

^ = -zi + g{{A + C)z, + {-A + C)z2) + g{A{z^ - z^)), (3) 

^ = -Z2 + 9{{~A + C)zi + {A + C)z2) + g(-A(zi - z^)) (4) 

and is a useful model system for gaining intuition about the behaviour of solu- 
tions of (IT]). The constant A is positive while C is non- negative. This system 
is symmetric under interchange of zi and Z2- Suppose that {x\,X2,x*^-,x\) is a 
stationary solution of (ll|) with the restrictions on the parameters atj leading to 
(|3|-([4|) and let z1 — x\ +a;5, z^ = x^ +x\. Then {z\^ Z2) is a stationary solution 
of ^-^ and 

xl=g{{A + C)zl + {-A + C)z*2), (5) 

xl=g{{-A + C)zl + {A + C)zl), (6) 

xl^g{A{zl-z^)), (7) 

xl^g{-A{zl-z*2)). (8) 

Conversely, if {zl.z^) is a stationary solution of p|-(|4]) and {x\,X2,x*^,x\) is 
defined by ©-(Is]) then a stationary solution of M is obtained. The fact that 
the starting solution was stationary assures the consistency conditions that z\ = 
x\ + X3 and Z2 = X2 -\- xX- Thus proving the existence of stationary solutions of 
([5|-([8l) also leads to a proof of the existence of corresponding stationary solutions 

3 Existence of stationary solutions 

In this section results are proved about the existence of stationary solutions of 
(JSI-Q. The techniques used to do so will be generalized in the next section to 
prove analogous results about stationary solutions of (fTl) . A picture which helps 
to provide an intuitive understanding of the results of this section and which 
played an important role in developing the theorems and their proofs is given 
as Fig. [11 The assumptions on the coefficients in (IT]) imply that C > in the 
system ^-(|4]) derived from it. Nevertheless the possibility that C is negative is 
allowed in this section. The reason for this is that dealing with negative values 




Figure 1: The nuUclines of ([Sj-Q in an example (A=3, C=1.65) 



of C in (|3|-(|4| yields intuition for the treatment of (II]) in the next section where 

it is allowed that some of the coefficients Cij are negative. For a positive real 

number e with e < 1 let Ki{e) = [0, e] x [2 - e,2], K2{e) = [2 - e, 2] x [0,e], 

K^{e) = [1 - e, 1 + e] X [2 - e, 2], K^ie) = [2 - e, 2] x [1 - e, 1 + e]. 

Theorem 1 (i) Suppose that ry+ < 1 is a positive constant and that -^ < ?7+. 

Then if A is sufficiently large there exists an e > such that there is a unique 

stationary solution of (Is]) and Q in each of Ki{t) and -f^2(e)- 

(ii) Suppose in addition that 77_ > ^ and ^ > ?;- . Then for A sufficiently large 

e can be chosen such that there is a unique stationary solution in each Ki{e) 

within 1,2,3,4. 

Proof Because of the symmetry of the system it is enough to prove the result 

for Ki and K^. Define a mapping by 



(zi,z2) ^{g{{A + C)zi + {- 
g{{-A + C)z, + {A- 



A + C)z2)+g{A{zi-z2)), 
^C)z2)+g{-A{z^-z2))). 



(9) 



Alternatively we can write (j){zi,Z2) — (0i(-Zi, ^2), </'2(^i, 22))- Stationary solu- 
tions of ([3]) and Q are in one to one correspondence with fixed points of the 
mapping (j). To prove the theorem it is enough to show that for A sufficiently 
large and a suitable choice of e the mapping leaves each Ki invariant and 
that the restriction of to each of these sets is a contraction. Suppose that 
e < j(l — ?7+). Let Aq be a positive constant and assume that A > Aq. The 
hypotheses of the theorem imply that —A + C < —(1 — rj^)A. Consider first 
the set Ki. For (zi, Z2) E Ki 

(A + C)zi + {-A + C)z2 < (1 + ii+)Ae - (1 - t^+)A{2 - e) 

= 2A[-(l-r,+ ) + e]<-(l-,7+)Ao. (10) 



Similarly A{zi — z^) < Ae— {2 — e)A < —Aq. Next the arguments of g occurring 
in the second component will be estimated. 

{-A + C)zi + {A + C)z2 > -A{1 + f]+)e + A{1 - r]+)(2 - e) 

= 2Ail-ri+-e)>{l-rj+)Ao (11) 

and —A(zi — Z2) > Aq. Note that 

l+tanhx= 7^, 1 — tanha;= t^. (12) 

Hence g{{A + C)zi + {-A + C)z2) < e-2(i-';+)^o and 5(^(21 - 23)) < er^^". 
If Aq is large enough then e can be chosen so that the inequality e^^^^^^+^^° < 
I is satisfied. Under these circumstances it follows that (j>i{zi,Z2), which is 
evidently positive, is less than e. In the same way gi{—A + C)zx + {A + C)z2) > 
1 - e-2(i-''+)^c) and g{-A{zi - 23)) > 1 - e'^^o. It follows that with the 
restrictions on Aq and e already assumed (1)2, which is evidently less than two, 
is no less than 2 — e. Hence Ki is invariant under 0. To prove that is a 
contraction it is helpful to use the identity g'{x) — 2g{x){\ ~ g{x)), which shows 
that g'{x) < 2min{e^^, e~^^}. It follows by the mean value theorem that if x 
and y both have the same sign and modulus greater than xq then 

\gix)-giy)\<2e-'-''\x-y\ (13) 

Hence if (2:1,22) and {■wi,W2) both belong to Ki then 

|(/.,(zi,Z2)-</),(z«i,ii;2)|<8(l + r7+)yle-2(i-')+)^(niax{|zi-Z2|,|«;i-W2|}) (14) 

for i = 1,2. The function xe^^f^^''+-'^ has its unique maximum on the interval 
[0,(^) at 2(1^ • Thus if Ao > 2(y^ and Aoe-2(i-')+)Ao < i(i + ^^)-i then 
the restriction of (/) to Ki is a contraction in the maximum norm. 

A similar analysis can be carried out in the case of K^ when the assumptions 
of part (ii) of the theorem hold. For (zi, Z2) G K^. 



{A + C)zi + {-A + C)z2 > -(1 + e)A + (3 - e)C 
' '3-e' 



>A{1 



1 



?]_ — 1 



(15) 



Suppose that e < §p^. Then it follows that {A + C)zi + {-A + C)z2 > 
''^~ Aq. On this region A{zi — Z2) < —\Aq. 

{-A + C)zi + (A + C)Z2 > -Ail + e)+ A{2 - e) > ^Aq (16) 

and -A{zi - Z2) > ^Aq. If it is assumed that e'^° < f and e~(3»j_-i)Ao ^ | 
these inequalities imply that Kg, is invariant. That the restriction of (f) to Kg is a 
contraction follows as in the case of Ki under the assumptions that ^0 ^ 3 — rj 
and Age^*^'^''-^^^"^" < |(1 + 77+)^^. This completes the proof of the theorem. 



The symmetry of the system imphes that the diagonal zi ~ Z2 is an in- 
variant submanifold. Stationary solutions on the diagonal are in one to one 
correspondence with points where ^{zi + Z2) = 0. When zi = Z2 

-(zi + Z2) - 1 - {zi + Z2) + 2g{C{zi + Z2)). (17) 

Thus stationary solutions on the diagonal are in one to one correspondence with 
zeroes of the function F{x) = 1 + 2g{Cx) — x on the interval (0,4). The facts 
that F(0) — 2 and F{A) < —1 show that F has a zero. The function x — 1 
is monotone increasing while for C < the function g{Cx) is monotone non- 
increasing. Thus in that case the zero of F is unique. For C > on the other 
hand, when the graph of 2g{Cx) crosses that oi x — \ its derivative must be 
no greater than one. Moreover in this case the second derivative of 2g{Cx) is 
negative and so the derivative of 2g{Cx) is less than one for all greater values of 
X. Thus its graph cannot cross that of cc — 1 again. It follows that the zero of F 
is unique in all cases and that there is always precisely one stationary solution 
on the diagonal. Hence in the cases where Theorem 1 shows the existence of 
four stationary solutions there are at least five in total. Note that if C is large 
and positive then the stationary point is close to (§,§)• 

Theorem 2 Under the hypotheses of Theorem 1 the system pl-Q has at least 
seven stationary solutions. 

Proof By Theorem 1 and the discussion following it there exist at least five 
stationary solutions. Using the symmetry of the system it is enough to show 
that there exists at least one other stationary solution. Consider the nuUclines of 
the system (|3])-(|4|, i.e. the zero sets of the right hand sides of the equations. On 
the first nullcline if zi is fixed the monotonicity of g implies that there is at most 
one possible value for Z2. Combining this with the implicit function theorem 
shows that the nullcline is the graph of a function 72 of zi defined on some open 
subset. On any closed subinterval of this subset the value of Z2 is bounded. It 
follows that in fact 72 is defined on the whole interval (0, 2). In a similar way it 
can be shown that the second nullcline can be written in the form zi — 71 (Z2) for 
a function 71 defined on (0, 2). The part of the curve zi — 71 (Z2) where zi < 1-l-e 
is contained in the strip defined by 2 — 2e~^^° < -22 < 2. This can be seen as 
follows. From estimates very similar to those derived in the proof of Theorem 
1 it can be concluded that for any fixed value of zi no greater that 1 -f e the 
mapping Z2 1— > ipzizi, Z2) defined on the interval [2 — e, 2] is a contraction. The 
unique fixed point of this mapping lies on the curve zi = 71(22) and satisfies the 
estimate claimed. Let zl j and Z2 i be the zi and Z2 coordinates of the stationary 
solution in Ki, respectively. To prove the theorem it is enough to show that in 
the interval {z^ 1, zj' 3) there is a value of Zi for which 72 is less than 2 — 26"^"^" 
and a value of zi for which 72 is greater than 2. For this it is helpful to compute 
the derivative of 72 . 

1 = [{-A + C)g'{{A + C)z, + {-A + C)z2) - Ag'{A{z^ ~ 7.2))]'^ 

dzi 

+ [{A + C)g'{{A + C)zi + {-A + C)z2) + Ag'{A{zi - Z2))] (18) 



The second expression in square brackets is small and, in particular, less than 
one half. Taking it onto the left hand side shows that -j^ is negative and that 

- ^ > ^ [(A - C)g\{A + C)zi + {-A + C)z2) + Ag' {,A{z^ - z^))]'^ (19) 



It can be concluded that 



dl2 



dzi 



> l^-le2(l-r,+ )Ao (20) 



in Ki and K3. This guarantees the existence of values of zi with the desired 
properties and completes the proof. 

When the linearization of the system about a stationary point in one of 
the Ki,l < i < A, is computed the contributions of the nonlinear terms are 
extremely small. As a consequence the linearization about a point of this type 
is very close to minus the identity. It follows that the stationary point is a 
hyperbolic sink. 

4 Stationary solutions of the full system 

In this section some of the techniques introduced in Sect. [3] are applied to 
obtain information about stationary solutions of (fTj) which covers an open set of 
parameters for this system whose definition is based on explicit inequalities. The 
system of equations defining stationary solutions can be simplified by replacing 
the variables Xi by diXi. This means that when proving theorems it can be 
assumed without loss of generality that di = 1 for all i. Corresponding results 
for general values of di can then be obtained by transforming the stationary 
solutions back to the original variables. First Theorem 1 will be generalized. 
Note that if a stationary solution of ^-^ is close to (0, 2), (2, 0), (1, 2) or (2, 1) 
then the corresponding stationary solution of (fTl) defined in Sect. [2] is close to 
(0,1,0,1)), (1,0,1,0), (1,1,0,1) or (1,1,1,0) respectively. For a positive real 
number e with e < 1 define 

Ki(e) = [0,e]x[l-e,l]x[0,e]x[l-e,l], 
^2(e) = [l-e,l]x[0,e]x[l-e,l]x[0,e], 
^3(e) = [l-e,l]x[l-e,l]x [0, e] x [l-e,l], 
Kiie) = [1 - e, 1] X [1 - e, 1] x [1 - e, 1] x [0, e]. 

Theorem 3 (i) Suppose that di = I for each i, that 77+ < 1 and A are pos- 
itive constants, that \cij\ < 77+|6ij| for all i and j and that maxij{|6ij|} < 
Aminij{|6ij|}. Then if minij{|6ij|} is sufficiently large there exists an e > 
such that there is a unique stationary solution of M in each of Ki{e) and K2{e). 
(ii) Suppose in addition that A— 1 is sufficiently small, 77- > g, Cij > r]- viisxk,i{\bki\} 
for i = 1,2 and all j. Then for nunij{|6y|} sufficiently large there exists an e 
such that there is a unique stationary solution in each Ki{e) with i = 1, 2, 3, 4. 



Proof The proof uses the same strategy as that of Theorem 1. Define a mapping 

: R* ^ R4 by 

(f)i{xi,X2, Xs, X4) == .9 X! "'J^J ^^^^ 

Stationary solutions of M are in one to one correspondence with fixed points of 
cj). To prove the theorem it suffices to show that each set Ki is invariant under 
(j) and that the restriction of (j) to each Ki is a contraction. In comparison with 
the proof of Theorem 1 the estimates here give less explicit information about 
the conditions in the statement of the theorem. Let Aq be a positive constant 
and assume that \bij\ > Aq for all i and j. Consider first the set Ki. Then 

J2 aijXj < -2(1 - 77+ - Me))Ao, (22) 

j 
J2 a2jXj > 2(1 - 7/+ - Me))Ao, (23) 

i 

^asjx, <-2(l-Afe))Ao, (24) 

j 

^ aijXj > 2(1 - Me))Ao (25) 

3 

for a constant M depending only on A, provided e is sufficiently small. These 
inequalities imply that Ki is invariant under (p if 

g2|e|g-4(l-r,+ -Me)Ao < ^ (26) 

and this is the case if ^o is sufficiently large and e is chosen appropriately. It 
can be shown that the restriction of to Ki is a contraction using a straight- 
forward generalization of the argument used in the proof of Theorem 1. Up to 
a numerical factor the contraction constant is 

e2|«lAmin{|6„-|}e-4(i-"+-*^^'™--{l^-l> (27) 

and this can be made as small as desired by choosing Aq large enough. The 
same arguments can be used to prove that K2 is invariant under ^ and that the 
restriction of 4> to K2 is a contraction. 

Next part (ii) of the theorem will be proved. This will be presented for 
the case of K^ - the case of K^ can be treated in the same way. Note that 
if ^ = 1 — A^^ then assuming that A — 1 is sufficiently small is equivalent 
to assuming that ^ is sufficiently small. Moreover for all i,j, fc, Z the quantity 



\btj\ - \bki\ 



is bounded by ^niaxi.j \bij\. In this case 
^ aijXj > (377_ - 1 - C - 6e)Ao, 

asjXj <^{l-2^-3e))Ao, 

Y^ a4.]Xj > (1 - 2^ - 3e))Ao. 



'J 



1,2, 



(28) 
(29) 
(30) 



This shows that for ^ and e sufficiently small the set K^ is invariant under (j). 
That the restriction to Ky, is a contraction can be shown as in the case of Ki . 

Note that for the parameter values which reduce (fTl) to (Isjl-Q the hypotheses 
of Theorem 3 reduce to those of Theorem 1. 

The argument used to obtain the stationary solution on the diagonal for 
the system pl-Q does not generalize to the system (IT]) since the latter has no 
symmetry. Instead the observation concerning the approximate position of that 
solution for A large made in the last section will be used. Define 



iiTs (6) = [l-e,l]x[l-e,l] 



1 


1 




'1 


1 





-e, - +e 


X 


— - 


-e, ^ +e 


2 


2 




2 


2 



(31) 



The direct analogue of the arguments used to treat Ki ioT 1 < i < A does not 
seem to work for K^ which is apparently not invariant under 0. To overcome 
this (j) will be replaced by another mapping if). 

Theorem 4 Assume that the hypotheses of Theorem 3 hold. Then if A — 1 is 
sufficiently small and miuij- \bij\ is sufficiently large there exists an e > such 



that there is a stationary solution of ([l| in Kr^. 
Proof Define new variables by ?/i = xi — 1, j/2 = 2:2 

J/4 = xa. 



\. Then 



1, ys 



X3, 



\ and 



bsixi + b'i2X2 + 6332^3 + bziXi 

= biiivi - y2) + ^33(2/3 - Vi) + {bsi + &32)2/2 + {b-ss 
1 



&34)2/4 



+ (&31 + &32)+ 2(^33 + ^34) 

641X1 + 6422:2 + 643X3 + 644X4 

= 641(2/1 - ^2) + 643(2/3 - 2/4) + (^41 + 642)2/2 + (643 + 644)2/4 

+ (&41+&42)+ 2(^43+^44) 

The first two terms in each of these expressions are obstructions to showing that 
(j) maps K^ into itself. To go further note that for any positive real numbers 
{x,w) 



g{x + w) 



1 



1/ 7 



(32) 



10 



This suggests rewriting the equation Xi = g{hi) for z = 3, 4 as 



l{K) 



1 



1. 



hiiiyi -2/2) - T^^iziy?, -Vi)- 



(33) 
(34) 



where 

1 

2 2"'"'"" "^' 2^ 

For i — 1,2 the equation is rewritten as y^ = Si where si = .9(^1) ^ 1- These 

equations are schematically of the form Ny — s for some matrix TV. Conditions 

will now be given which ensure that the matrix N is invertible. For this it is 

"J " 

^21 ^22] 
N22 is invertible then N is also invertible and the inverse of N is given by 

Now consider the matrix N22- Its determinant is 



useful to write it as a matrix 



of two-by-two blocks. If the matrix 



'^2V^^21 



N. 



22 



1, 



1-2^33 



1, 



^43 



1 



^33^43 



-"33 



^33 



^33 



2^3-3') 



1 - 2633I - 



"33 



^43 



1 - 26, 



33 , 



- 1 



-4(1 - h 



33 I 



2 

^33 



^33 
633 + ^43 



- 1- 



"33 



^43 



"33 



"43 



"33 



(35) 



Let /3 = mini J l^ijl- Then 633^ < /3 ^. It follows in particular that if /3 > 2 then 



1-6 



33 



> 



In addition 



&33+fa43 



< ^. It follows that the square bracket in 

(35) is bounded above by —2(1 — ^) and from below by —4(1 -I- £). In can be 
concluded that the following inequality holds 



/3A(l+0<detiV22< 



1 



(1-0- 



(36) 



In particular, the determinant does not vanish and so N is invertible. Moreover 
the modulus of the inverse of the determinant can be bounded by a constant 
multiple of /3~^. It follows that the norm of N22 ^ can be bounded by a constant 
and hence that the norm of N~^ can be bounded by a constant multiple of /3. On 
the translate ofSs by (— 1,— 1,— i,— i) define a mapping by ^{y) = N^^s{y). 
For i = 1,2 



M> (2r;_-2^-5e)(min|6,, 



3,4 



(37) 



Choose ^ < |?7_ and e < jq'H-- Then the first bracket on the right hand side 
of this equation is greater than 2^-- This in turn gives exponential bounds for 
Si and 32- It remains to estimate S3 and S4. Only the argument for S2 will be 
given since the corresponding argument for S4 is very similar. 



1^31(2/1 - 2/2) + ^33(^3 -yi)\ < 3/3Ae. 



(38) 
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Using this and ( 33 1 gives the estimate 



|s3|<i[9/32A2e2+4/3e] 



(39) 



By choosing ^ and f3'^e sufficiently small it can be ensured that S3 and S4 are 
bounded by an arbitarily small multiple of e. Now e can be chosen to satisfy 
this condition and to ensure at the same time that si and S2 are bounded by 
an arbitrarily small multiple of e. In this way it can be ensured that ip leaves 
i?5 invariant. The conclusion of the theorem follows by applying the Brouwer 
fixed point theorem to ^. 

Remark This theorem also holds, with the same proof, if the hypothesis on 77^ 
is weakened to 77- > 0. 

It has proved possible to extend the method used to obtain Theorem 4 to 
obtain the existence of analogues of the remaining stationary solutions of p|- 
Q. Unfortunately this has only been achieved for a quite restricted choice of 
parameters. To understand where the condition for the parameters comes from, 
note that if 



Ci 

C2 



-612 



C12 



Ci4 



611+Cii 

-hi - C21 - 624 - C24 



^22 



C22 



then 



5(/ii(Ci,l,0,l))=.g(/i2(l,C2,l,0)) = - 



(40) 
(41) 

(42) 



The cases treated in the next theorem are those where Ci or (,2 is close to ^. 
Let 

X [l-e,l]x[0,e]x[l-e,l], 
-e,-+e X [l-e,l] X [0,e]. 
Theorem 5 Assume that the hypotheses of Theorem 4 hold and that 

Cii + 2ci2 + 2ci4 



^6(e) = 


"1 1 

2"^'2 + 


e 


kjie) = 


l-e,l]x 


'1 
2 



/3 



-3 



<e 



(43) 



If ^ is small enough and /3 is sufficiently large then e can be chosen so that there 
is a stationary solution of ([l| in Ke{e). An analogous statement holds if the 
inequality ( 43 1 is replaced by 



C21 + 2C22 + 2C24 



^ 



<e 



(44) 



and K^le) is replaced by K'j(e) in the conclusion. 
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Proof The strategy of the proof follows that of Theorem 4. It is enough to 
treat the case of KQ{e) since that of Ki{e) is very similar. This time define new 
variables by yi = xi - 5, 2/2 = 2^2 - 1, 2/3 = 2:3 and yi = x^- 1. 

(611 + cii)xi + (612 + Ci2)x2 + (613 + Ci3)a;3 + (614 + Ci4)x4 

= (foil + Cii)yi + (612 + Ci2)y2 + (&13 + Ci3)y3 + (614 + Ci4)y4 

-611 +612 + foi4+ 2^11 +C12 +ci4J (45) 

In this case the equation xi — g{hi) will be rewritten as 

yi-^{bii+cii)yi--{bi2 + ci2)y2-^{bi3+ci'i)y^--{biA + ci4)yi = si (46) 
where 

Sl = 9{hl) - 2 ~ 2 *-^" ^"^"^^^ - 2 (^12 +Cl2)y2 - 2 (^13 + 213)2/3 - 2 (^14 + Cl4)2/4 

(47) 
For i = 2,3,4 the equation Xi = g{hi) is written as y2 = 5(^12) — 1, 2/3 = .9(2/3) 
and 2/4 — gihi) — 1. So as in the proof of Theorem 4 the condition for stationary 
solutions has been rewritten in the form Ny = s for a matrix N. In this case 

where iVn is scalar with the only entry 1 — 2011, 



N is of the form 



/ 

A'^12 is a row matrix with entries —^au, i = 2,3,4, and the identity matrix is 
three by three. 1 — ^on < 1 — -^(3. If /3 > 3 the matrix is invertible. The inverse 
of iVii can be bounded by a constant multiple of /3~^. The inverse matrix N~^ 
is easily computed. Its norm can be bounded by a constant multiple of /3. On 
the translate of Bq by (— ^, —1, 0, — l) define a mapping by ^jJ{y) — N^^s{y). It 
will now be shown that ^ maps its domain into itself. The estimates 

h2 > (1 - e - 3e)/3, (48) 

hs <-(l-^-3e)/3, (49) 

hi > (1 - C - 3e)/? (50) 

suffice to take care of the last three components Sj . It remains to treat the first 
component. Once this has been done the estimate for N^^ which is already 
available completes the proof. The quantity si can be estimated using the 



inequality ( 32 1 



I (fell + cii)2/i + (612 + 012)2/2 + {bi3 + ci3)2/3 + (^14 + ci4)yi| < 8/3Ae. (51) 

It follows that 

|si|<32/32aV+Qa + i')/3C (52) 

This estimate leads to the desired conclusion as in the proof of Theorem 4. 
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For the special choices of the parameters leading to the system (Isll-Q the 
quantities Ci and C2 reduce to 



In this case, as ^ varies from one to ^, the Q vary from zero to one. Each 



of the extra conditions (|43| and (|44|) on the parameters reduces to ^ = |. 



When Theorem 5 applies this gives some information about the position of the 
stationary solutions whose existence is guaranteed by Theorem 2. 



5 Conclusions and outlook 

The dynamical system introduced in i4] to describe the interactions between T 
cells and macrophages has been analysed with respect to its stationary solutions. 
For certain open sets of the parameter space it was shown that there are at least 
seven stationary solutions of which at least four are stable. The positions of 
these stationary solutions were described approximately. The same features are 
already found in the model system (Isll-Q with the two parameters A and C. 
There it is only necessary to assume that ^ < ^ < 1 and that A is sufhcently 
large to get these conclusions. Under weaker assumptions, which reduce to 
^ < 1 and A sufhciently large for the model system it was shown that there 
are at least three stationary solutions of which at least two are stable. By 
contrast it was shown in [7] that under a suitable smallness assumption on the 
parameters there is a unique stationary solution and that all solutions converge 
to it at late times. For the model system a sufficient condition implying the 
latter behaviour is 2^1 + |C| < 1. 

As already remarked in [7], for C < A the system pl-Q is competitive 
and hence any solution converges to a stationary solution as a consequence 
of the results of [2]. In particular there can be no periodic solutions. This 
conclusion depends on the fact that the system is two-dimensional and so cannot 
be drawn for the system (II]) , even in the case of parameter values for which it is 
competitive. The results of this paper do not rule out the possibility that there 
are values of the parameters for which the model system has periodic solutions. 
They also do not prove that there cannot be more than seven stationary solutions 
or more than four stable stationary solutions for some values of the parameters. 
Simulations have given no indications that either of these phenomena (periodic 
solutions or extra stationary solutions) actually occur. 

From the biological point of view one of the most interesting new findings 
of this paper is the coexistence of stationary solutions such as those close to 
(0, 1, 0, 1) and (1, 1, 1, 0) with the same values of the parameters. Certain small 
changes of the parameters can move the system from a situation where both of 
these exist to a situation where one of them has disappeared. On a heuristic 
level this gives a scenario where a small external influence on the system leads 
to a large change in its properties. Note that these two stationary solutions 
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have very different biological properties. In the second case the Thl and Th2 
cells are secreting a comparable large amount of cytokines. The total cytokine 
concentration is 50% greater in the second case than in the first. The model of 
[3] is crude in may ways but their are many models for biological systems which 
share qualitative features with this one and the kind of switching phenomenon 
which has been observed here may occur much more widely. Thus it is good to 
have mathematical tools which allow it to be studied. 

The Thl/Th2 picture of autoimmune diseases is no longer up to date. Newer 
models involve players other than Thl and Th2 cells. In particular Thl7 cells 
have come to play a very important role [S]. Mathematical insights obtained 
for one model can be useful in understanding other models, even those which 
remain to be invented. Thus the fact that the biological content of the model 
of [4] is not close to the explanations of the behaviour of the immune system 
which are presently most popular does not prevent it being worthwhile to study 
its mathematical properties in detail. 

The mathematical description of the immune system in the model of [5j is at 
the level of the interaction of populations of cells and what happens in a given 
cell is handled as a simple black box. Another very interesting task is to model 
the processes within cells which determine whether they differentiate into types 
such as Thl and Th2. These usually involve the behaviour of transcription 
factors. Some models of this kind are studied in [TT] and [T]. Models of the 
population of cells and models of the molecular machinery within cells describe 
systems whose intrinsic nature are very different. Nevertheless they may be re- 
lated mathematically and exploiting this kind of relation is a particular strength 
of mathematics. 
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